podatki <- read.table("/cloud/project/Poglavje 5/Naloga 1/Medicinsko osebje.csv", header=TRUE, sep=";", dec=",")
head(podatki)
##   BolnišnicaID OddelekID OsebaID Stres Starost Izkušnje Spol Oddelek
## 1            1        11       1     7      36       11    0       0
## 2            1        11       2     7      45       20    0       0
## 3            1        11       3     7      32        7    0       0
## 4            1        11       4     6      57       25    1       0
## 5            1        11       5     6      46       22    1       0
## 6            1        11       6     6      60       22    1       0
##   BolnišnicaVelikost Usposabljanje
## 1                  2             1
## 2                  2             1
## 3                  2             1
## 4                  2             1
## 5                  2             1
## 6                  2             1

Opis spremenljivk:

podatki$SpolF <- factor(podatki$Spol, 
                        levels = c(0, 1), 
                        labels = c("Moški", "Ženska"))

podatki$OddelekF <- factor(podatki$Oddelek, 
                           levels = c(0, 1), 
                           labels = c("Splošni", "Spec"))

podatki$VelikostF <- factor(podatki$BolnišnicaVelikost, 
                            levels = c(1, 0, 2),
                            labels = c("Srednja", "Majhna", "Velika"))

podatki$UsposF <- factor(podatki$Usposabljanje, 
                         levels = c(0, 1), 
                         labels = c("NE", "DA"))

str(podatki)
## 'data.frame':    1000 obs. of  14 variables:
##  $ BolnišnicaID      : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ OddelekID         : int  11 11 11 11 11 11 11 11 11 12 ...
##  $ OsebaID           : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Stres             : int  7 7 7 6 6 6 6 7 7 6 ...
##  $ Starost           : int  36 45 32 57 46 60 23 32 60 45 ...
##  $ Izkušnje          : int  11 20 7 25 22 22 13 13 17 21 ...
##  $ Spol              : int  0 0 0 1 1 1 1 1 0 0 ...
##  $ Oddelek           : int  0 0 0 0 0 0 0 0 0 1 ...
##  $ BolnišnicaVelikost: int  2 2 2 2 2 2 2 2 2 2 ...
##  $ Usposabljanje     : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ SpolF             : Factor w/ 2 levels "Moški","Ženska": 1 1 1 2 2 2 2 2 1 1 ...
##  $ OddelekF          : Factor w/ 2 levels "Splošni","Spec": 1 1 1 1 1 1 1 1 1 2 ...
##  $ VelikostF         : Factor w/ 3 levels "Srednja","Majhna",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ UsposF            : Factor w/ 2 levels "NE","DA": 2 2 2 2 2 2 2 2 2 2 ...
fit_ML <- glm(Stres ~ Starost + Izkušnje + SpolF + OddelekF + VelikostF + UsposF, 
              data = podatki)

library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:DescTools':
## 
##     Recode
vif(fit_ML)
##               GVIF Df GVIF^(1/(2*Df))
## Starost   2.987514  1        1.728443
## Izkušnje  2.990715  1        1.729368
## SpolF     1.003130  1        1.001564
## OddelekF  1.001023  1        1.000511
## VelikostF 1.001278  2        1.000319
## UsposF    1.000940  1        1.000470
podatki$StdOst <- round(rstandard(fit_ML), 3)
podatki$CooksD <- round(cooks.distance(fit_ML), 3)

hist(podatki$StdOst,
     main = "Histogram standardiziranih ostankov",
     ylab = "Frekvenca",
     xlab = "Standardizirani ostanki")

podatki$StdOcene <- scale(fit_ML$fitted.values)

library(ggplot2)
ggplot(podatki, aes(y=StdOst, x=StdOcene)) +
  theme_linedraw() +
  geom_point() +
  xlab("Standardizirane ocenjene vrednosti") +
  ylab("Standardizirani ostanki")

head(podatki[order(podatki$StdOst), c("OsebaID", "StdOst")])
##     OsebaID StdOst
## 304     304 -3.700
## 306     306 -2.793
## 301     301 -2.693
## 27       27 -2.578
## 824     824 -2.562
## 823     823 -2.552
head(podatki[order(-podatki$StdOst), c("OsebaID", "StdOst")])
##     OsebaID StdOst
## 107     107  2.578
## 8         8  2.498
## 635     635  2.413
## 804     804  2.399
## 896     896  2.358
## 652     652  2.219
head(podatki[order(-podatki$CooksD), c("OsebaID", "CooksD")])
##     OsebaID CooksD
## 28       28  0.012
## 304     304  0.011
## 27       27  0.010
## 896     896  0.010
## 21       21  0.009
## 22       22  0.009
podatkiC <- podatki[c(-304),  ]

fit_ML <- glm(Stres ~ Starost + SpolF + OddelekF + VelikostF + UsposF, 
              data = podatkiC)

summary(fit_ML)
## 
## Call:
## glm(formula = Stres ~ Starost + SpolF + OddelekF + VelikostF + 
##     UsposF, data = podatkiC)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.19070  -0.54182  -0.02076   0.56696   2.19272  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      5.956317   0.115439  51.597  < 2e-16 ***
## Starost         -0.004565   0.002143  -2.130   0.0334 *  
## SpolFŽenska     -0.449683   0.058427  -7.697 3.37e-14 ***
## OddelekFSpec     0.079331   0.051542   1.539   0.1241    
## VelikostFMajhna -0.489366   0.056332  -8.687  < 2e-16 ***
## VelikostFVelika  0.430683   0.076298   5.645 2.16e-08 ***
## UsposFDA        -0.735682   0.051566 -14.267  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.6634201)
## 
##     Null deviance: 942.64  on 998  degrees of freedom
## Residual deviance: 658.11  on 992  degrees of freedom
## AIC: 2434.1
## 
## Number of Fisher Scoring iterations: 2
library(Rmisc)
## Loading required package: lattice
## Loading required package: plyr
## -------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## -------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename,
##     summarise, summarize
## The following object is masked from 'package:ggpubr':
## 
##     mutate
Statistika <- summarySE(podatkiC, 
                        measurevar="Stres", 
                        groupvars=c("BolnišnicaID"), 
                        conf.interval=0.95)

library(ggplot2)
ggplot(Statistika, aes(x=BolnišnicaID, y=Stres )) +
 geom_bar(position=position_dodge(), 
           stat="identity") +
 geom_errorbar(aes(ymin=Stres-ci, ymax=Stres+ci),
               width=0.1,                    
               position=position_dodge(.9)) +
 theme_linedraw()

Statistika <- summarySE(podatkiC, 
                        measurevar="Stres", 
                        groupvars=c("OddelekID"), 
                        conf.interval=0.95)

ggplot(Statistika, aes(x=OddelekID, y=Stres )) +
 geom_bar(position=position_dodge(), 
           stat="identity") +
 geom_errorbar(aes(ymin=Stres-ci, ymax=Stres+ci),
               width=0.1,                    
               position=position_dodge(.9)) +
 theme_linedraw()

library(lme4)
library(lmerTest)
hierar_fit0 <- lmer(Stres ~ (1 | BolnišnicaID) + (1 | OddelekID), 
                 REML = FALSE,  
                 data = podatkiC)

summary(hierar_fit0)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: Stres ~ (1 | BolnišnicaID) + (1 | OddelekID)
##    Data: podatkiC
## 
##      AIC      BIC   logLik deviance df.resid 
##   1941.3   1960.9   -966.6   1933.3      995 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.98593 -0.66090  0.01731  0.67055  3.12364 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  OddelekID    (Intercept) 0.4840   0.6957  
##  BolnišnicaID (Intercept) 0.1588   0.3985  
##  Residual                 0.2992   0.5470  
## Number of obs: 999, groups:  OddelekID, 100; BolnišnicaID, 25
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   5.0026     0.1072 24.9966   46.66   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ICC_Bolnica <- 0.159/(0.484+0.159+0.299)
ICC_Bolnica
## [1] 0.1687898
ICC_Oddelek <- 0.484/(0.484+0.159+0.299)
ICC_Oddelek
## [1] 0.5138004
hierar_fit1 <- lmer(Stres ~ Starost + SpolF + OddelekF + VelikostF + (1 | BolnišnicaID) + (1 | OddelekID), 
                 REML = FALSE,  
                 data = podatkiC)

summary(hierar_fit1)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: 
## Stres ~ Starost + SpolF + OddelekF + VelikostF + (1 | BolnišnicaID) +  
##     (1 | OddelekID)
##    Data: podatkiC
## 
##      AIC      BIC   logLik deviance df.resid 
##   1815.8   1859.9   -898.9   1797.8      990 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1493 -0.6789  0.0123  0.6770  3.4268 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  OddelekID    (Intercept) 0.4801   0.6929  
##  BolnišnicaID (Intercept) 0.0537   0.2317  
##  Residual                 0.2610   0.5109  
## Number of obs: 999, groups:  OddelekID, 100; BolnišnicaID, 25
## 
## Fixed effects:
##                   Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)       5.506546   0.156513  63.725971  35.183   <2e-16
## Starost          -0.002432   0.001392 905.793627  -1.747   0.0810
## SpolFŽenska      -0.438058   0.038346 907.598306 -11.424   <2e-16
## OddelekFSpec      0.060560   0.142344  75.153566   0.425   0.6717
## VelikostFMajhna  -0.486770   0.187205  24.937236  -2.600   0.0154
## VelikostFVelika   0.418833   0.245375  25.046486   1.707   0.1002
##                    
## (Intercept)     ***
## Starost         .  
## SpolFŽenska     ***
## OddelekFSpec       
## VelikostFMajhna *  
## VelikostFVelika    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
## Warning in abbreviate(rn, minlength = 6): abbreviate used with non-
## ASCII chars
##             (Intr) Starst SplFŽn OddlFS VlksFM
## Starost     -0.383                            
## SpolFŽenska -0.187  0.016                     
## OddelekFSpc -0.454 -0.001 -0.002              
## VelkstFMjhn -0.512 -0.002  0.002  0.000       
## VelikstFVlk -0.392 -0.001  0.006  0.000  0.327
hierar_fit2 <- lmer(Stres ~ Starost + SpolF + OddelekF + VelikostF + UsposF + (1 | BolnišnicaID) + (1 | OddelekID), 
                 REML = FALSE,  
                 data = podatkiC)

summary(hierar_fit2)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: 
## Stres ~ Starost + SpolF + OddelekF + VelikostF + UsposF + (1 |  
##     BolnišnicaID) + (1 | OddelekID)
##    Data: podatkiC
## 
##      AIC      BIC   logLik deviance df.resid 
##   1788.8   1837.9   -884.4   1768.8      989 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1321 -0.6946 -0.0210  0.6790  3.3633 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  OddelekID    (Intercept) 0.31674  0.5628  
##  BolnišnicaID (Intercept) 0.09455  0.3075  
##  Residual                 0.26108  0.5110  
## Number of obs: 999, groups:  OddelekID, 100; BolnišnicaID, 25
## 
## Fixed effects:
##                   Estimate Std. Error         df t value
## (Intercept)       5.859025   0.162311  71.478220  36.097
## Starost          -0.002529   0.001391 907.790344  -1.818
## SpolFŽenska      -0.437505   0.038316 910.365779 -11.418
## OddelekFSpec      0.061044   0.117169  74.675107   0.521
## VelikostFMajhna  -0.486579   0.187213  24.937276  -2.599
## VelikostFVelika   0.419107   0.245386  25.046570   1.708
## UsposFDA         -0.697899   0.117174  74.688560  -5.956
##                     Pr(>|t|)    
## (Intercept)          < 2e-16 ***
## Starost               0.0694 .  
## SpolFŽenska          < 2e-16 ***
## OddelekFSpec          0.6039    
## VelikostFMajhna       0.0155 *  
## VelikostFVelika       0.1000    
## UsposFDA        0.0000000786 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
## Warning in abbreviate(rn, minlength = 6): abbreviate used with non-
## ASCII chars
##             (Intr) Starst SplFŽn OddlFS VlksFM VlksFV
## Starost     -0.373                                   
## SpolFŽenska -0.179  0.016                            
## OddelekFSpc -0.360 -0.001 -0.003                     
## VelkstFMjhn -0.494 -0.002  0.002  0.000              
## VelikstFVlk -0.378 -0.001  0.006  0.000  0.327       
## UsposFDA    -0.364  0.009 -0.004  0.000  0.000  0.000
hierar_fit3 <- lmer(Stres ~ Starost + SpolF + OddelekF + VelikostF + UsposF + UsposF:SpolF + (1 | BolnišnicaID) + (1 | OddelekID), 
                 REML = FALSE,  
                 data = podatkiC)

summary(hierar_fit3)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: 
## Stres ~ Starost + SpolF + OddelekF + VelikostF + UsposF + UsposF:SpolF +  
##     (1 | BolnišnicaID) + (1 | OddelekID)
##    Data: podatkiC
## 
##      AIC      BIC   logLik deviance df.resid 
##   1790.5   1844.5   -884.2   1768.5      988 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1278 -0.6972 -0.0085  0.6743  3.3587 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  OddelekID    (Intercept) 0.31619  0.5623  
##  BolnišnicaID (Intercept) 0.09464  0.3076  
##  Residual                 0.26102  0.5109  
## Number of obs: 999, groups:  OddelekID, 100; BolnišnicaID, 25
## 
## Fixed effects:
##                        Estimate Std. Error         df t value
## (Intercept)            5.844102   0.164351  75.122236  35.559
## Starost               -0.002543   0.001391 907.833863  -1.828
## SpolFŽenska           -0.415969   0.053725 909.340597  -7.743
## OddelekFSpec           0.060847   0.117074  74.655910   0.520
## VelikostFMajhna       -0.486721   0.187190  24.937258  -2.600
## VelikostFVelika        0.419264   0.245356  25.046525   1.709
## UsposFDA              -0.665757   0.129878 112.325951  -5.126
## SpolFŽenska:UsposFDA  -0.043820   0.076638 910.408046  -0.572
##                      Pr(>|t|)    
## (Intercept)           < 2e-16 ***
## Starost                0.0678 .  
## SpolFŽenska          2.59e-14 ***
## OddelekFSpec           0.6048    
## VelikostFMajhna        0.0154 *  
## VelikostFVelika        0.0999 .  
## UsposFDA             1.24e-06 ***
## SpolFŽenska:UsposFDA   0.5676    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Warning in abbreviate(rn, minlength = 11): abbreviate used with non-
## ASCII chars
## 
## Correlation of Fixed Effects:
## Warning in abbreviate(rn, minlength = 6): abbreviate used with non-
## ASCII chars
##             (Intr) Starst SplFŽn OddlFS VlksFM VlksFV UspFDA
## Starost     -0.365                                          
## SpolFŽenska -0.237 -0.001                                   
## OddelekFSpc -0.355 -0.001 -0.004                            
## VelkstFMjhn -0.488 -0.002  0.001  0.000                     
## VelikstFVlk -0.373 -0.001  0.005  0.000  0.327              
## UsposFDA    -0.392  0.000  0.301 -0.001 -0.001  0.000       
## SplFŽn:UFDA  0.159  0.018 -0.701  0.003  0.001 -0.001 -0.433
anova(hierar_fit0, hierar_fit1, hierar_fit2, hierar_fit3, test="Chi")
## Data: podatkiC
## Models:
## hierar_fit0: Stres ~ (1 | BolnišnicaID) + (1 | OddelekID)
## hierar_fit1: Stres ~ Starost + SpolF + OddelekF + VelikostF + (1 | BolnišnicaID) + (1 | OddelekID)
## hierar_fit2: Stres ~ Starost + SpolF + OddelekF + VelikostF + UsposF + (1 | BolnišnicaID) + (1 | OddelekID)
## hierar_fit3: Stres ~ Starost + SpolF + OddelekF + VelikostF + UsposF + UsposF:SpolF + (1 | BolnišnicaID) + (1 | OddelekID)
##             npar    AIC    BIC  logLik deviance    Chisq Df
## hierar_fit0    4 1941.3 1960.9 -966.64   1933.3            
## hierar_fit1    9 1815.8 1859.9 -898.89   1797.8 135.5013  5
## hierar_fit2   10 1788.8 1837.9 -884.41   1768.8  28.9599  1
## hierar_fit3   11 1790.5 1844.5 -884.25   1768.5   0.3268  1
##                Pr(>Chisq)    
## hierar_fit0                  
## hierar_fit1     < 2.2e-16 ***
## hierar_fit2 0.00000007389 ***
## hierar_fit3        0.5675    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1